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Abstract 

A  random  walk  in  a  discrete  space  of  large  dimensions  is 
devised  to  find  the  free  energy  of  finite  Ising  spin  systems  with 
nearest  neighbor  interactions.   The  system  studied  is  the  three 
dimensional  cubic  lattice,  finite  and  periodic  in  two  of  the 
dimensions,  and  infinite  in  the  third  dimension.   The  effect  of 
low-order  surface  terms  in  the  graphical  expansion  of  the  free 
energy  is  shown  by  removing  them  from  the  numerical  results  for 
the  finite  lattice,  and  comparing  this  reduced  free  energy  to  the 
bulk  free  energy.   Using  transfer  matrix  formalism,  a  Monte  Carlo 
solution  is  obtained  by  establishing  a  correspondence  between  the 
eigenvalue  equation  satisfied  by  the  transfer  matrix,  and  the 
collision  density  equation  of  the  random  walk.   The  walk  is  then 
simulated  to  provide  Monte  Carlo  estimates  of  the  free  energy, 
magnetization,  and  spin  correlation. 


I.   Introduction 

One  important  area  of  investigation  in  physics  at  the  present 
time  is  that  of  cooperative  phenomena,  where  correlations  between 
particles  are  so  strong  that  standard  independent  particle  approxi- 
mations fail  to  give  accurate  predictions  of  system  behavior. 
Unfortunately,  even  oversimplified  models  have  proven  extremely 
difficult  to  solve.   This  thesis  will  look  at  one  of  these  models, 
the  Ising  model,  in  which  spin  one-half  particles  are  fixed  .to 
lattice  sites,  and  will  investigate  how  the  size  of  the  lattice 
affects  one  of  its  properties,  the  Helmholtz  free  energy.   This  will 
be  done  by  comparing  numerical  results  for  various  size  lattices 
with  the  predictions  of  graphical  expansions. 

A  new  Monte  Carlo  method  which  is  based  on  the  transfer  matrix 
technique  is  developed  to  numerically  solve  this  model.   By 
sampling  a  random  walk  having  the  proper  statistical  features, 
estimates  of  the  free  energy,  magnetization,  and  spin  correlation 
are  obtained. 

A.  The  Ising  Model 

Consider  a  system  of  particles,  with  spin  one-half,  fixed  at 
the  vertices  of  a  lattice.   At  each  site,  associate  a  binary  scalar 
quantity,  a,  which  can  be  thought  of  as  the  projection  of  the  spin. 
The  spin  can  be  thought  of  as  pointing  in  some  predefined  direction 
in  space  (say,  the  +z  direction)  if  a   =  +1  and  in  the  opposite 
direction  if  o  =  -1. 

Define  the  interaction  energy  between  particles  as  follows : 


(1.1)     E. .  =  -a. .Jo   o      ;   a..  =  1  if  spin  i  and  spin  j  are 

nearest  neighbors 
'a.  .  =0  otherwise 


That  is,  if  two  nearest  neighbors  are  pointed  in  the  same 
direction,  E. .  =  -J;  and  if  they  are  pointed  in  opposite  directions 
E   =  +J.   Notice  that  J  is  a  measure  of  the  strength  of  the 
coupling,  and  J  >  0  favors  parallel  spins  by  giving  them  a  lower 
interaction  energy. 

Consider  also  that  the  system  may  interact  with  an  external 

magnetic  field.   If  each  spin  has  a  magnetic  moment  y  ,  its 

s 

interaction  with  a  field,  H  =  H£  %,    will  be 


(1.2)  (Eu) .  =  -y  Ha. 

H'j     "s   j 


The  full  Hamiltonian  for  a  lattice  of  N  spins  will  be 


,   N    N  N 

(1.3)     "HT  =  "  k  5~'  J~~  a.  .Jo.o  .    -   y  H  J~~  o  . 


or 


N 
(1.4)    T4-  =  -  ZZ  Jo±o.    -   ysH  YZ  o. 

<ij>    J       j=l   J 


wher«»  "Nij>  means  sum  over  all  nearest  neighbor  pairs. 

Using  the  above  Ising  Hamiltonian,  one  can  immediately 
write  the  canonical  partition  function 


(1.5) 


z  =  Z 


B 


where  the  sum  is  over  all  system  configurations,  a  configuration 
being  a  particular  set  of  the  values  of  all  the  spins.   T  is  the 
absolute  temperature;  kg  is  the  Boltzmann  constant.   Equation  (1.5) 
can  also  be  expressed  as: 


(1.6) 


Z  -  I 


zz    zz  ••• 

^_          exp 

> 

a1=±l   o2=±l 

°N=±1 

<ij 

Ka.a.+  h>   a  . 
.   1  J    ^1-3 


(1.7) 


where  K  =  6J  ;   h  =  By  H 
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Notice  that  this  problem  can  be  reformulated  in  the  language 
of  quantum  mechanics.   Interpret  the  a's  as  the  Pauli  matrices  a 
To  be  precise 


(1.8) 


(1.9) 


©  ©  ..<j) ...  $FT)  © 


a.  =1  x  1  x...a  ...x  1  x  l 


1   0 


* 


1  = 


1   0 


0   1 


) 


where  direct  product  notation  has  been  used.   That  is,  choose  the 
direct  product  of  all  the  spins  of  the  lattice,  with  the  a  's  dia- 
gonal, as  the  representation.   A  general  basis  set  will  be 
written  as 


(1.10) 


|u>.  | 


a1>o2> 


°N>; 


=  ±i 


and  for  a  particular  ket  the  following  notation  will  be  used 

(1.11)  !++++...,> 

The  quantum  Hamiltonian  has  the  exact  form  as  (1.5), 
since  it  is  diagonal  in  this  representation,  but  the  a's  must 
be  interpreted  as  operators.   In  this  representation  equation 
(1.5)  becomes  identically: 


-B'W- 
(1.12)  Z  =  Tr[e    I] 


The  equilibrium  thermodynamic  variables  follow  in 
standard  fashion.   The  free  energy  per  spin: 


k  T 
(1.13)  F  =  -  -jj-  An  Z 


The  internal  energy  per  spin: 

2 
(1  iin  U  =  -  ^-  3(£n  z) 

and  the  magnetization  per  spin: 

N   ^ 

(1.15)       m  =  N°l  ¥      =  1  »(*"  Z) 

N         N    3h 


B.  Physical  Reality  of  the  Model1 
Mathematically,  the  problem  is  one  of  associating  with  each 
lattice  site  a  binary  variable,  which  for  convenience  has  been 
taken  as  ±1.   A  Hamiltonian  has  been  introduced  which,  as  will 
be  seen,  favors  ordered  states. 

This  model  can  represent,  in  a  simple  way,  the  following 
three  physical  systems. 

1.   A  magnetic  system,  where  each  binary  variable  re- 
presents the  projection  of  the  spin  angular  momen- 
tum of  a  particle  fixed  at  the  site. 


(1.16)  sz  =  | 


At  the  time  of  the  model's  introduction,  it  was 
hoped  that  the  Ising  Hamiltonian  would  lead  to  an 
explanation  of  spontaneous  magnetization,  where  all 
or  most  of  the  spins  align.   The  language  of  spin 
systems  was  used  in  Section  IA  to  introduce  the  model 
and  will  be  used  throughout  this  thesis. 
A  mixture  of  two  kinds  of  molecules,  say  A  and  B. 
a  =  +1  would  mean  that  the  site  is  occupied  by  an  A 
atom,  and  a  =  -1  would  mean  that  the  site  is  occupied 
by  a  B  atom.   This  interpretation  could  lead  to  a 
description  of  ordering  in  binary  alloys;  for  example, 
for  J  >  0,  the  Hamiltonian  favors  A  atoms  grouping 
with  A  atoms,  and  B  atoms  grouping  with  B  atoms. 


3.   A  mixture  of  atoms  and  holes,  for  which  a  =  +1 
means  that  the  site  is  occupied  by  an  atom,  and 
a   =  -1  means  that  the  site  is  vacant.   Interpreted 
in  this  manner,  the  model  is  called  a  lattice  gas, 
and  could  be  used  as  a  model  for  condensation. 

The  correspondence  between  the  thermodynamic  variables  of 
the  above  three  systems  can  be  made  rigorous.   In  particular, 

the  isotherm  1/l(h)  for  the  magnetic  problem  can  be  easily  trans- 

2 
formed   into  P(v),  the  pressure  versus  the  specific  volume  for 

the  lattice  gas  problem. 

The  approximate  aspects  of  the  model  are  the  following: 

1.  fully  localized  spins  (or  atoms) 

2.  nearest  neighbor  interactions  only 

3.  high  anisotropy  (spin  system  only). 
However,  cooperative  phenomena  are  thought  to  be  relatively 
independent  of  the  explicit  form  of  the  Hamiltonian,  especially 
at  the  critical  point.    Also,  the  Ising  Hamiltonian  must  be 
understood  before  proceeding  to  more  physical,  and  correspondingly, 
more  difficult  problems. 

C.   Order-Disorder  Interpretation  of  the  Model 

To  show  how  ordered  states  are  preferred,  let  us  look  at 

the  model  from  a  statistical  mechanical  viewpoint.   At  zero  field, 

the  probability  of  the  system  being  in  the  state 

(1.17)  |u>  =  |a1a2, .  .  •crN)>   ,   a  =  ±1 


IS 


(1.18) 


Pr{ |u>  }  =  Z 


-1 


+J  XZ  a   a    /kRT 


'B' 


At  high  temperatures,  the  denominator  of  the  exponent  dominates 
All  the  states  are  nearly  equiprobable ,  becoming  so  as  T  ■+  °°. 
As  the  temperature  is  lowered,  states  which  are  highly  ordered 


(i.e.,  large  value  of  J  ^Z  o.a.)    have  a  higher  probability 

<ij>  X    J 
than  those  that  are  less  ordered.   The  lower  the  temperature, 

the  greater  the  difference  in  the  probability  o.f  two  different 

states.   Notice  that  this  type  of  order  does  not  prefer  either 

the  plus  or  minus  z  direction. 

Upon  adding  a  magnetic  field,  the  probability  becomes 


(1.19)     Pr{ |u>  }  =  Z 


-1 


+J   JViW 
<1J>   J 


+  uH  XZ  ai/kBT 
J 


,1 


+  uH  ^T"  o./kRT 

1~  J 
The  new  factor  [e  ]   breaks  the  symmetry,  making 

more  probable  those  states  with  spins  pointing  predominantly  in 

the  +z  direction. 

The  above  picture  of  a  system  becoming  more  ordered  as  the 

temperature,  is  lowered  does  not,  however,  suggest  either  the 

presence  or  the  absence  of  a  phase  transition. 


D.  A  Brief  History  of  the  Model 
In  view  of  the  voluminous  literature  on  the  Ising  model, 
only  the  high  points  in  its  development  will  be  covered.   The 


n 

model  was  first  formulated  by  Ernst  Ising   in  1925.   In  his 
dissertation,  suggested  by  his  mentor  Wilhelm  Lenz,  Ising  exactly 
solved  the  model  in  one  dimension,  finding  no  phase  transition. 

The  next  important  advance  came  in  1941  when  Kramers  and 
Wannier   showed  that  the  Ising  partition  function  was  the  largest 
eigenvalue  of  a  matrix,  called  the  transfer  matrix.   By  means  of 
a  symmetry  argument,  they  found  the  critical  point  (K  =.44069) 

of  the  2D,  oo  x  oo}  square  lattice. 

c 
The  above  lattice  with  h  =  0  was  fully  solved  by  Onsager 

in  1942,  who  followed  and  extended  the  method  of  Kramers  and 

Wannier.   He  obtained  the  same  critical  temperature,  a  continuous 

energy,  and  a  specific  heat  that  became  infinite  as  -log|T-T  |, 

where  T   is  the  critical  temperature. 

A  derivation  of  the  spontaneous  magnetization  for  the 

7  ft 

lattice  was  published  in  1952  by  C.N.  Yang  .   Griffiths  , 

in  1968,  showed  that  Yang's  result  is  a  lower  bound  to  the  true 

magnetization.   An  important  set  of  inequalities  was  also 

o 
developed  by  Griffiths  .   One  of  them  states  that  for  the  zero 

field  case  the  covariance  of  any  pair  of  spins  is  nonnegative 
(J>0);  that  is,  the  probability  that  they  have  the  same  sign  is 
at  least  one-half.   The  above  lattice  in  a  finite  magnetic  field 
has  yet  to  be  exactly  solved;  in  3D  nothing  is  known  rigorously. 
Many  approximate  methods  have  been  used  on  the  Ising  model, 
in  particular,  the  self-consistent  field  approach   .   Since  the 
cooperative  phenomena  are  true  many-body  problems,  these  approxi- 
mations did  not  reproduce  the  correct  behavior  of  the  thermo- 
dynamic variables.   However,  the  use  of  exact  series  expansions 


8 


has  proved  extremely  fruitful.   A  correspondence  can  be 
established  between  the  exact  power  series  expansion  of  the 


partition  function  and  combinatorial  problems  in  graph  theory 
For  example,  Wakefield    obtained  the  first  twelve  terms  of 
the  high  temperature  series  for  the  zero  field,  simple  cubic 


12 
lattice.   Knowledge  of  these  terms  has  led  to  good  approximations 

of  the  thermodynamic  variables. 


II.   Formulation  of  the  Problem 

A.  Transfer  Matrix  Formalism 
Basic  to  most  problems  in  equilibrium  statistical  mechanics 
Is  the  partition  function  Z,  from  which  the  free  energy  per 
particle  (spin)  can  be  obtained: 


k  T 
(2.1)  F  =  -  -f-  In    " 


It  will  now  be  shown  that  the  Ising  partition  function  can  be 

written  as  the  largest  eigenvalue  of  a  matrix  called  the  transfer 

13 
matrix.   Among  the  many  proofs  available,  that  of  Kac    will  be 

followed,  for  it  is  both  the  simplest  and  the  most  elegant. 

First,  the  concept  of  layer  must  be  defined.   A  lattice  can 
be  thought  of  as  being  composed  of  layers  of  sublattices.   For 
example,  in  1-D  a  single  spin  is  a  layer.   For  the  2-D  square 
lattice,  having  n  rows  and  m  columns,  a  row  of  spins  is  a  layer. 
Likewise,  the  3-D  simple  cubic  lattice  of  size  m  *  n  *  £,  has 
a  2-D  square  lattice  (mxn)  as  a  layer. 

There  are  two  obvious,  but  important  facts.   First,  the 
Ising  Hamiltonian  couples  nearest  neighbor  layers.   Second, 
knowing  the  configuration  of  all  the  layers  uniquely  fixes  the 
configuration  of  the  entire  lattice. 

The  following  derivation  will  be  made  explicitly  using  the 

square  lattice.   The  argument  extends  trivially  to  the  simple 

cubic  lattice  with  only  a  redefinition  of  layer  and  changes  in 

summation  variables. 

Denote  the  configuration  of  the  jth  layer,  u..   For  example, 

J 

10 
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for  the  square  lattice  with  n  rows  and  m  columns:, 

(2.2)         u.  =  {a      ar),.  .  .o    }      ;       a   =  ±1 

the  ordered  set  of  all  spin  directions  in  the  jth  row.   The 
partition  function  becomes 


(2.3) 


■$\ 

=  >         a 

all  system 
configurations 

A 

■  > 

all-laver  1 
conf igs . 

( conf igs . ) 

> 

) 

> 

all  layer  2  all  layer j  '"'all  layer  n 
conf igs.     conf igs.        conf igs 


(2.4)       Z  =  H   E!  !  •••  511  +iL 


-BX-j 


ul    u2 


Create  periodic  boundary  conditions  by  wrapping  the  lattice 
on  a  cylinder  so  that  the  spins  of  then' throw  interact  with  the 
spins  of  the  first  row,  and  then  bend  the  cylinder  into  a  torus 
by  having  the  spins  of  the  first  column  interact  with  the  spins 
of  the  last  column.   Then 


(2.5)          Z  =  >  V    V     .  .  .V  „ 

v   J'              L u-,u0  uou0    u  u, 

u, jUp,...u  12   23     n- 1 

where  the  2   dimensional  matrix  V  is  defined  as: 


m  m  m 

(2.6)    <u|V|u'>    =   V  =    exp[K  ^Z  a,a,+1+h   YZ  a,+K  YH  o.a\l 

3=1      J    J  j=l      J         ,1=1      J    J 
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where  the  first  term  in  the  exponential  is  the  intralayer 
contribution,  the  second  is  the  interlayer  contribution,  and  the 
last  contribution  is  from  the  field,  where 


(2.7)  am+1  =  o±    , 


and  where  a  direct  product  basis  of  all  the  spins  in  a  row  is  used: 

|u>  =  |a1,a2,o3,...am> 

(2.8)  |u'>  =  |a^,a',c',...a^> 
o,c'    =  ±1  . 

Since  all  indices  are  dummy  indices 

(2.9)  Z  =  YZ   (Vin.  )n  -  Tr(Vn)  =  H"  X? 

u  i 


where  the  A.'s  are  the  eigenvalues  of  V 
i  & 


(2.10)  V|TX>  =    Ai|<F1>   ;       \1    >    X2    >      X.. 


Now   let   the   number  of  rows    in   the    lattice   become   infinite 


(2.ii)  z  =  YZ  a?  =  x?[i+  Or1  )nJ 

i        x  X  i>l    Al 


(2.12)  Lim    (^—-)    =    *nX, 

n->oo  n  -1 


12 


Ik 
Furthermore,  it  can  be  shown    that  the  probability  that  the  row 

j  is  in  the  state  ju.^is 


(2.13)     pH|u>  }  =  Ku.lv1>!2 

J  J 


B.   Factorization  of  the  Transfer  Matrix 
It  has  been  shown  that  the  partition  function  per  row  for 
an  m  by  infinite  lattice  is  the  largest  eigenvalue  of  the  matrix 

V,  the  transfer  matrix.   V  can  be  put  into  a  more  convenient  form 

15 
by  factoring  it  into  the  product  of  simpler  matrices   .   Define 

the  following  2   dimensional  matrices 


m 

C2.14)        (V0)     =  <5   ,  exp[K  >  'cc..,] 
•  2  uu'     uu1    ^   4 — =-  j  j+1 

J=l   J  d 

m 

(2.15)  (V_)   ,  =  6   ,  exp[h  J~'   a.] 

v  3  UU'      UU'     K     Vry   J 

m 

(2.16)  (V,,n'  E   exP^K  ^Z  a  a'.] 

-LUU  ,•  _  t    J   J 


consider 


(2'17)        Wl]uu'  =  ^Vuu"<V3W^lW 

u  u  ' 


(2.18)  =  (V0)   (V-)   (V,  )   , 

v     '  v  2'uuv  3  uuv  1  uu' 
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(2.19) 


m_  m_  m 

exp[K  YZ  ff-iai+i+h   ZZ  CTi+K  5      a  . a  '.  ] 


J=l 


j=l      d  j=l 


(2.20) 


V       , 


Therefore 


(2.21) 


v  =  v2v3Vl 


V      can   be    put    in   the    following    fo 


rm 


(2.22)  (Vuu«    =    exP[K   XZ  tf-cr'.]    =   TT'exp[K   o  .a}^ 


J    J  J     J 


J  J 


Looking  at  one  factor  in  the  (direct)  product 


(2.23) 

(2.24) 
where 


exp[K  a .o \ ] 
J  J 


'  K    -K  \ 
e     e    \ 


-K 
e     e 


K  J 


_  -,  K     x  -K 

=  le   +  a  e 


(2.25)         1 


/  1   0  \ 


0   1 


0   1  \ 


a  o 


Therefore 


(2.26) 


V1  =  7J  [e  +a  e   ] 

j 


where 


111 


Y     ©    ©  . . .  Ql  .  •  .©>  © 

(2.27)        a.  =  1  x  1  x...a  ...*  T  x  1 

J 


Rewriting 


(2.21)  V  =  V2V3V1 


where 


(2.28)        V2  =  TT  exp[K  oj.aj  +  1] 

J 


(2.29)  V^  =  TT  exp[ha  ] 

5  J        J 

(2.30)  V  =  TT  [eK+a*e"K] 

1    J       J 


Notice  that 


(2.3D       [V2,V3]  =  0  ;    CV^]  *   0  ;   [V^]  *  0 


and  that  a   is  a  spin  flipping  operator 


(=•32)         (;)-xC);    a)=°x'vo) 


C.   The  Character  of  the  Largest  Eigenvalue  and  Its  Corres- 
ponding Eigenfunction 
One  general  remark  should  be  made  about  the  2   dimensional 
transfer  matrix.   A  theorem  of  Frobenius    states  that  for  a 
finite  dimensional  matrix  with  all  elements  non-negative,  the 
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largest  eigenvalue  is  non-degenerate,  and  the  corresponding 
("largest")  eigenvector  has  non-negative  elements.   For  finite 
m,  this  theorem  is  obeyed  by  the  transfer  matrix,  since  all 
its  elements  are  of  exponential  nature. 

D.  High  Temperature  Series  Expansions 
1.   Introduction 
High  temperature  series  expansions  for  tne  Ising  free 
energy  will  be  derived  for  the  bulk  lattice  and  the  finite 
lattice.   It  will  be  seen  that  the  lowest  order  differences 
between  them  are  the  results  of  "surface"  terms  which  contribute 

to  the  finite  lattice,  but  not  to  the  bulk  lattice. 

17 
Following  Fisher   ,  the  procedure  for  the  expansion  will 

be  reviewed.      On  setting  the  magnetic  field  equal  to  zero, 

equation  (1.12)  becomes 

(2.33)         Z  =  Tr[TT  exp(K  o.a  )] 

where  the  product  is  over  all  nearest  neighbor  bonds.   Expanding. 
the  exponential  in  a  power  series,  one  sees  that 

(2.3*0        exp[K  a. a. J  =  (cosh  K)(l+o.c.\;)  ;   v  =  tanh  K 

so  that  equation  (2.33)  becomes 


(2.35)         Z  =  coshqN/2K  Tr[JT   (1+o.a.v)] 

<ij>      J 
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where  N  is  the  number  of  spins;  qN/2  is  the  number  of  bonds;  q 
is  four  for  the  square  lattice  and  six  for  the  simple  cubic 
lattice . 

Now  consider  the  expansion  of  the  product.   The  coefficients 
of  v   will  be  all  combinations  of  all  the  different  nearest 
neighbor  bonds  taken  £  at  a  time. 
However,  for  any  i 


(2.36)  Tr[(c.)P]  =  0    when  P  is  odd 

=  2    when  P  is  even. 


Therefore,  all  terms  with  an  odd  number  of  a.'s  will  vanish 
under  the  trace.   The  surviving  terms  can  be  put  in  one-to-one 
correspondence  with  graphs  drawn  on  the  lattice.   These  graphs 
will  be  closed  polygons  since 

a)  from  the  nature  of  the  expansion  of  the  product,  each 
bond  can  occur  only  once  in  each  term 

b )  to  survive  the  trace  an  even  number  of  bonds  must  meet 
at  any  lattice  point. 

Therefore  one  can  write 

00 

(2.37)  Z    =    2N(cosh  QN/2K)[1+   Yl'   PpOOv    ] 

£  =  0 


where  Pn(N)  is  the  number  of  different  graphs  having  ft  bonds 
that  can  be  drawn  on  a  lattice  of  N  sites. 

2.   The  Bulk  Lattice 
The  first  few  P  ' s  of  the  square  bulk  lattice  can  be  found 

by  inspection. 
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(2.38)         P1(N)  =  P2(N)  =  P3(N)  =  0 


PI|(N)  =  N>    P5(N0  =.  0,     Pg(N)  =  2N  . 

Notice  that  Pn(N)  would  have  contributions  both  from  octagons 
and  separated  squares. 

By  formally  taking  the  logarithm,  one  can  show  that  the 
series  expansions  for  the  free  energy  per  spin  become 


(2.39)      -  -r-^m  =  in  2  +  |  £n  cosh  K  +  Yl     pov£ 

where  P.  is  the  coefficient  of  the  term  in  Pp (N)  which  is  linear 
in  N.   That  is 


(2.40)         P£(N)  =  NPJ  +  N2P^  +, 


This  is  to  be  expected  from  the  intensive  nature  of  the  bulk 
free  energy  per  spin. 

3.  The  Finite  Lattice 
Consider  a  finite,  periodic  square  lattice,  of  width  m, 
wrapped  on  a  cylinder  of  infinite  length.   All  graphs  of  less 
than  m  bonds  would  contribute  in  the  same  way  to  both  the  finite 

-l  o 

and  bulk  lattices.   However,  when  I    equals  m,  surface  graphs 
appear  on  the  finite  lattice  but  not  on  the  bulk  lattice.   They 
are  rings  which  circumscribe  the  cylinder,  and  they  contribute 


18 


terms  of  order  v   to  the  expansion  of  the  partition  function 
of  the  finite  lattice.   Higher  order  surface  terms  (of  order 
vm   ,  vm   etc)  can  also  be  found. 
Therefore,  to  leading  order  in  v 


F      F 
(2.41)  b  _    f  =  c(m)vm 

KB      B 


or 


F      F 
(2.42)         FT  "  ITT  =  C(m)exp[-m|£n  v|] 
B      B 


where  F,  is  the  free  energy  per  spin  of  the  bulk  lattice,  Ff 
is  the  free  energy  per  spin  of  the  finite  lattice,  and  C(m)  is 
an  unknown  function  of  m. 

Equation  (2.42)  shows,  to  lowest  order,  the  size  dependence 
of  the  free  energy  per  spin. 

To  see  this  in  more  detail,  contrast  the  expansion  of  the 
bulk  lattice  with  the  finite,  3  x  ro  ,  periodic  lattice. 


(2.43)         Zb(N)  =  2N  coshqN/2K[l  +  P^(N)vlj+P6(N)v6+...] 


(2.44)         Zf(N)  =  2N  coshqN/2K[l+S3(N)v3  +  Pi|(N)v1|+S5(N)v5 

6.o  ^  7+...] 


+Pg(N)v  +S?(N)v 


where  S„(N)  is  the  number  of  surface  graphs  having  I   bonds 
that  can  be  drawn  on  a  lattice  of  N  sites.  - 
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Therefore 


i. 


q 


(2.45)  -  ! — fft  -  tn  2  -  4  In  cosh  K 
kBT  2 


P2  8 
lrD  4X1D  6XD  8  ^4V   , 
N[V  +P6V  +P3V  "  -2—  + 


and 


(2.46) 


*B 


T 


-    £n    2    -    %   in    cosh    K    = 


i[S3v3  +  Pi)vll  +  S5v5  +  P6v6  +  S7v7 


c2    6 
2~ 


s3y 


7 


+  ...] 


For  this  lattice  then 


S-.V- 


(2.47) 


b     f   _   3 


kBT    kBT 


+  0(v5) 


For  the  periodic  simple  cubic  lattice  of  size  m  x  m  x  «>,  the 
low  order  surface  terms  were  found  to  be 


(2.48) 


S  (N) 
m 


=  2/m 


WN) 


=  4(m-l) 


(2.49) 


S^(N)   =  I2(m-1)  +  Mm-l)(m-2)(m-3)  +  12(m_1)(m_2) 
N  3 
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III.   The  General  Theory  of  Monte  Carlo 

19 
As  a  formal  definition  of  Monte  Carlo,  that  of  Shreider 

will  be  used: 


The  Monte  Carlo  method  is  a  method  of  solving 
various  problems  in  computational  mathematics  by 
constructing  for  each  problem  a  random  process  with 
parameters  equal  to  the  required  quantities  of  that 
problem.   The  unknowns  are  determined  approximately 
by  carrying  out  observations  on  the  random  process 
and  by  computing  its  statistical  characteristics 
which  are  approximately  equal  to  the  required  para- 
meters . 


A  simple  example,  which  has  every  feature  of  the  above 
definition,  is  the  integral 

b 

(3-D  I  =  J  L(x)dx 

a 


Divide  and  multiply  the  integrand  by: 


(3-2)  AZJ.  where  P(x)  >  0;  c  =     P(x)dx 


then 


(3.3)  I  -  ej  $}  <^>  a, 


If  {P(x)/c}  is  interpreted  as  a  probability  density,  then 


(3.4)  I  -   «<$f}> 
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that  is,  the  integral  I  is  c  times  the  average  value  of  the 

L(x ) 
function  p-? — y  ,  with  respect  to  the  density  P(x)/c. 

Suppose  that  one  can  obtain  a  sequence  of  numbers 

20 


following  function,  !„„„,  is  an  unbiased  stastistical  estimate 


x-,,Xp...xN,  each  sampled  from  the  density  P(x).   Then    the 

is  i 

L(x) 
of  the  mean  value  of  C^-t — ( 

P(x) 

N   L(x  )/P(x 

(3.5)         i  :  c  n  — h E  W  =  c  < 


N         -  XEST     \P75FT, 


The  samples  x, ,  x?...xN  need  not  be  statistically  independent 
This  particular  form  of  estimator  is  called  the  sample  mean. 
Observe  the  trivial  case  where  P(x)  =  L(x): 


L(x.)/L(x.)       N 


(3.6)     IEST  -  C  }_  -        N        -  c  N  -  c  , 

i  =  l 


making  the  estimator  exact  (zero  variance).   As  suggested  by 
the  above  trivial  case,  the  statistical  error  6 


(3.7)         6  =  |I-IEST 


depends  on  the  choice  of  P(x).   Great  increases  in  accuracy 
(efficiency)  can  be  obtained  by  choosing  P(x)  such  that  the 
ratio  p4 — y  remains  as  constant  as  possible  over  the  range 
a,b.   This  procedure  is  called  importance  sampling.   From  the 
law  of  large  numbers,  the  error  also  depends  on  N,  the  size  of 
the  sample.   Other  important  theorems    (the  central  limit 
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theorem  and  the  Chebyshev  inequality)  can  be  invoked  to  show 
that,  when  the  samples  are  independent,  the  variance  of  the 
estimator  is  given  by 


(3.8) 


v(IEST)a  |  , 


and  therefore  the  standard  deviation  of  I 


EST' 


(3.9) 


S.D.  =  /v 


S.D.   a 


/N 


By  increasing  the  sample  size,  N,  one  can  make  the  standard 
deviation  as  small  as  one  likes.   As  one  can  estimate   the  mean 
value  of  a  function,  I,  one  can  also  estimate  the  standard 
deviation  of  I.   In  this  thesis  the  estimated  standard  will  be 
taken  as  the  measure  of  error. 

An  estimator  such  as  the  sample  mean,  equation  (3-5),  is 
a  random  variable  because  it  is  a  function  of  random  numbers.   It 
therefore  has  a  density  function  associated  with  it.   If,  as  is 
the  situation  with  the  sample  mean,  the  mean  of  this  new  random 
variable  equals  the  desired  answer,  the  estimator  is  called  an 
unbiased  estimator.   The  difference  between  the  mean  of  the  above 
density  function  and  the  desired  answer  is  called  the  bias. 


IV.   The  Monte  Carlo  Method 

A.  The  Relationship  Between  Random  Walks 

22-23 
and  Eigenvalue  Equations 

The  largest  eigenvalue  of  the  transfer  matrix  and,  there- 
fore,the  Ising  partition  function,  will  be  obtained  by  simulating 
a  random  walk  (a  Markov  chain  with  stationary  transition 
probabilities)  which  has  the  proper  statistical  features.   The 
walk   will  be  conducted  in  a  discrete  space  having  the  same 
dimensionality  as  the  transfer  matrix.   A  descriptive  language 
will  be  used  where  points  will  be  moved  from  one  "position"  in 
this  space  to  another. 

Consider  the  following  eigenvalue  equation  with  the  operator 
A,  the  eigenvalue  A  ,  and  the  eigenfunction  |  <f>> 


(4.D  I4>>  =  f|<j)> 


In  a  discrete  space  A  is  a  matrix;  in  a  continuous  space  A  is 
an  integral  operator.   In  what  follows,  the  notation  of  the 
discrete  case  will  be  used. 

Consider  first  the  special  case  where  the  operator  (-r-)  is 
stochastic;  that  is, 


(4.2)         ^±Hi  >  o,  all  u,u'  ;  }      (^1)  =  lf    an  u» 

X  u    \ 


(Notice  that  the  discrete  space  is  labeled  using  the  same  indices 
as  the  direct  product  spin  space  used  for  the  transfer  matrix.) 
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Equation  (4.1)  may  now  be  given  a  probabilistic  interpretation. 
Let  \$>,    the  eigenfunction ,  be  a  probability  density  function: 


(4.3)         Pr{u}  =  <u|<|)>  =  <f) 


'i 


That  is,  the  probability  of  the  point  being  at  position  u  is  tj>  . 
Because  equation  (4.1)  is  linear  and  homogeneous,  the  probability 
can  be  taken  as  normalized: 


(4.4)         7~  <u|(j»  =  1  . 
u 

A 
Also,  let  (r-)   f  be  the  conditional  probability  of  u,  given  u' 

'  A     UU  ' 


(4.5)  Pr{u|u'}    =    <u|    £    |u'>    =    (£)uu, 


That  is,  the  probability  of  the  point  moving  to  u,  given  that  it 

•  A., 
'  uu 


is  at  u',  is  (t-)„„i  •   Notice  that  equation  (4.4)  gives  the 


required  normalization.   Equation  (4.1)  can  be  written  as 


(4.6)      4>  =  IZ(x)^'<i>u, 

u    u,  A      u 


and  interpreted  as  follows: 

The  probability  of  being  at  u  after  the  move,  <j>  , 
equals  the  probability  of  being  at  u'  before  the  move, 
fa    ',  times  the  probability  of  moving  to  u  from  u',  (x")uu»> 
summed  over  all  u'. 

h    is  therefore  the  stable,  steady  state  density  for  a  point,  or 
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a  population  of  points  which  moves  by  means  of  the  "kernel", 


A 
X' 


This  process  can  easily  be  simulated.   Given  a  population 


of  points,  with  positions  u',  distributed  according  to  <j>  ,  ,  move 
each  of  them  independently  to  a  new  position,  u,  by  sampling 
(y)   , .   The  process  of  all  points  in  the  population  making  a 
move  will  define  one  iteration.   Ey  equation  (4.6),  the  population 
after  the  iteration  will  still  be  distributed  according  to  <£>.   How- 
ever, since  both  ^  and  X  are  not  known,  a  procedure  of  convergence 
must  Dt-  found.   More  about  this  important  consideration  will  follow 
The  eigenvalue  equation  can  be  generalized  to 


(4.7) 


> 


where 


A  is  stochastic,  Y-  A   , 

'  c —  uu' 


1;  A   ,  >  0  all  u,u' 
'   uu'  — 


B  is  positive  and  diagonal,   B   , 

^  b     '    uu' 


6   ,B  :  B   ,  >  0   all  u,u' 
uu'  u    uu '  — 


and  A  is  the  eigenvalue  of  the  operator  BA 
Rewriting 


B 


(4.8) 


u 


A   uu' 


u 


The  new  matrix,  r-,  can  be  incorporated  into  the  above  scheme  by 

A 

introducing  the  following  fissioning  technique.   Create  an 
arbitrary  set  of  probability  density  functions  g.,(0  on   the  space 
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of  non-negative  integers  £  such  that 


<"-9)       <$\u)    "  V*  E  Fu 


This  mean  value  will  be  called  the  fissioning  part  of  the  matrix. 
Use  them  in  the  following  way.   When  a  point  moves  from  u'  to  u, 
sample  g  U)  for  a  K    .   Fission  the  point  such  that  it  becomes  £ 
points  at  u.   Then  move  all  these  "daughters"  independently.   For 
a  proof  that  this  procedure  will  simulate  equation  (4.7),  see 
Appendix  A. 

One  last  generalization  of  the  eigenvalue  equation  is  needed 

Consider : 


(4.10)  |4>>  =  3£|<|>>   , 


where  now, 

B  is  still  diagonal  and  positive 

A  is  positive  (not  necessarily  stochastic) 

Rewriting, 


B  A   . 
u  uu' 


u' 


Choose  an  arbitrary  conditional  probability,  puu.>  whose  norm 
is  known: 


(4.12)        Puu,  >  0  all  u,u<  ;     YZ   Puu-  =  V 
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that  is, 


(4.13)  p^[u|u'}  =  P   ,/C  , 

v  '        uu'   u ' 


Multiply  and  divide  equation  (4.11)  by  this  conditional  probability 


(4.14)     *    .^^.^     ^i>, 

ll'  llll'    I      ll' 


B  C  ,  A  ,    r  P   , 
uu  '  J   ui1 ' 

uu '  L  u 


Now  consider  the  random  walk  which  arises  from  this  equation. 

P   , 
Move  the  points  from  u'  to  u  by  sampling" 


u ' 


Upon  arrival  of  a  point  at  u,  fission  the  point  to  £ 
daughters  such  that 


B  C  ,  A   , 
u  u'   uu' 


4.15)        <£>  =  ~^~-  t^-   =  P   , 

A    P.     uu' 


UU' 


The  right-hand  side  of  this  equation  will  again  be  called  the 

fissioning  part  of  the  matrix. 

Because  of  the  freedom  one  has  in  choosing  P   , ,  many 

b  uu"     J 

different  random  walks  have  the  same  density  |cj>>.   This  is  just 
importance  sampling  (see  Chapter  III)  applied  to  linear  equations, 
and,  here  also,  great  increases  in  accuracy  and  stability  can  be 
had  by  choosing  P   ,  such  that  the  fissioning  part  of  the  matrix, 

B  C  ,  A   , 
(4.16)  u  u   uu 


X        P   ,   » 
uu' 


remains  as  constant  as  possible 
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B.  Convergence 
The  relationship  between  stable  random  walks  and  eigen- 
value problems  was  established  in  the  previous  section.   Since 
neither  the  eigenfunction  nor  the  eigenvalue  is  known,  a 
method  is  needed  to  provide  convergence  (relaxation)  from  a 
trial  eigenfunction  to  the  true  eigenfunction.   Let  <f>.,  be  the 
vector  particle  density  after  the  Jl'th  iteration.   <}>   then  is  the 
trial  eigenfunction,  here  regarded  as  a  density  function. 
Also,  let  A„  be  the  trial  eigenvalue.   Consider  the  following  game 

1.  Sample  a  population  of  points  from  <p,,;  %   =    0. 

2.  Move  and  fission  each  point  in  the  population 
independently,  as  described  in  the  previous 
section.   From  the  discussion  following  equation 
(4.6),  it  can  be  seen  that. 


<*■«'  ^+i  =  -x^ 


That  is,  if  the  population  was  distributed 
according  to  p?  before  the  move,  it  is  dis- 
tributed  according  to  <J>n+-i  after  the  move. 
3.  I  ^>  I   +  1;  go  to  step  2  until  I   >_   ^MAX;  then 

go  to  step  4. 
k.      Stop. 
Looking  at  the  analytic  features  of  the  game,  let  <j)   and  A.  be 
the  eigenvectors  and  eigenvalues  of  BA  ordered  such  that 
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(4.18)        xi  >  A2  >   V  v 


Expand  c}>   in  eigenfunctions  of  BA : 


(4.19)        ?0  =  ZZ   C^1   ,       C.  =  (b1^ 

i 


Then,  from  equation  (4.17), 


A.   ^  . 

(4.20)  4>     =  TZ  Mt^H1 

1  T 

and 

(4.21)  %l   =   IZc/^V 


Let    Xm   =    A .  .       Then 

T  l 


"i      x —  ,xi,a„  ?i 


(4.22)  <L    =    C-.^    +   ^_    (t±)    C J> 

*  1  =  2       Al 


Assuming  no  degeneracy, 


(4.23)  Lim  J£  =  C^1 


Therefore,  for  large  I    (after  relaxation),  the  above 
game  produces  a  poplulation  of  points  drawn  from  <f)  ,  the 
"largest"  eigenfunction  of  BA . 

The  speed  of  convergence  can  be  seen  to  depend  on: 

1.  The  size  of  C, ,  i.e.  how  good  the  trial  eigenvector  is 

2.  The  ratio  of  \      to  A-,. 
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Those  readers  familiar  with  standard  numerical  methods  can 
see  that  the  above  method  is  the  stochastic  analog  of  the 
Neuman  iteration  method. 

Essential  to  the  above  proof  are  the  conditions 

(4.24)         <u|(j)1>   10,     all  |u>  , 

as  these  numbers  are  interpreted  as  probabilities.  The 
"largest"  eigenfunction  of  the  transfer  matrix  has  been 
shown  in  Section  II. C   to  satisfy  these  conditions. 

C.  Applying  Monte  Carlo  to  the  Transfer  Matrix 
The  concepts  of  Section  IV.  A  will  be  applied  to  the 
eigenvalue  equation 

(2.10)  |Y>  =  £|<F>  , 

where  V  is  the  transfer  matrix  and  X    is  its  largest  eigenvalue. 
By  suitably  factoring,  V  may  be  put  in  the  form  of  BA 
(equation  (4.7)  or  (4.10)).   The  discrete  space  in  question  is 
the  direct  product  spin  space  of  all  the  spins  in  a  layer  (row, 
for  2-D).   Points  then  will  move  from  layer  configuration  |u'^ 
to  layer  configuration  |u).   This  is  conceptually  equivalent 
to  having  spins  in  |u')>  flip  so  that  the  configuration  becomes 
|uy  .   From  the  explicit  form  of  the  transfer  matrix  it  will 
be  seen  that  this  is  indeed  what  happens. 
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1.  The  Obvious  Walk 

It  has  been  seen  in  Section  IV.  A  that  different  random 
walks  have  the  same  density  |<t^>  .   Also,  the  method  of  creating 
and  simulating  them  was  shown.   Because  of  the  previous  splitting 
of  V  into  three  matrices,  VpV_V, ,  VpV_  being  diagonal,  an  obvious 
random  walk  can  be  constructed  using  V,  as  the  moving  factor  and 
VpV_  as  the  fissioning  factor.   Because  of  its  pedagogical  value, 
and  because  it  was  the  first  walk  computed  for  this  thesis,  it 
will  be  examined  in  detail. 

The  transfer  matrix  for  a  2-D  square  lattice  with  m  spins 
per  row  is 


(2.21)  V  =  V2V3V1  , 


where 


(2.28)  V2  =  ]Texp(Kajaj+1) 

J 


(2.29)  V   =  TTexp(ha  )  , 

J       J 


and 


(2.30)  V     =  TT[eK+a^e'K] 

j      J 


K.  -K, 


After  dividing  and  multiplying  V,  by   T(e  +e  )  and  definin 


(4.25)  P  =  (e2K+l)  1   ,     0  <  P  ■:  1   , 
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V  takes  the  form 


(4.26) 


ni 


V  =  (2  cosh  k)m  TT[(l-P)+a*P] 

J=l       J 


Now  remove  the  number  (2  cosh  k)   from  V   and  place  it  with 
V„  so  that  new  operators  are  defined  as 


(4 


m-i 


V2  =  (2  cosh  k)  TJexp(Ka  a    ) 


m 
27)    <    Vn  =  JT  C(l-P)+a^P] 

I        .1=1        J 


S 


v  =  v2v3Vl 


J 


Since  V  V  is  diagonal,  it  has  the  properties  of  the  B  matrix 
of  the  second  walk  in  Section  IV.  A: 


(4.28) 


V2V3  <=£  B 


Now  the  matrix  V-,  must  be  checked  to  see  if  it  is  stochastic; 
if  so,  it  corresponds  to  the  matrix  A  of  the  previous  section 
V  is  stochastic  if 


(4.29) 


y  /u|V  |u'/  =  1    for  any  |u'^> 

u 


Consider  the  expansion  of  the  product 


in 


(4.30) 


V   =  IT  C(l-P)+a*P] 
j  =  l        J 
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Remembering  that  ax  is  a  flip  flop  operator,  that  the  | u V s 

are  orthogonal,  and  using  combinatorial  theory,  the  coefficient, 

a,  of  the  |  uj>  'th  component  of  V  |  u  'J>  is 


(4.3D  a  -  (l-P)m"NFD  PNFD   , 


where  NFD  is  the  number  of  spin  flips  necessary  to  change  u1 
into  u.   Therefore, 


(4.32)       ^IVju^  =  (1-P)m-NFD  PNFD 


Now  the  number  of  kets  which  differ  from  |u'J>  by  NFD  spins  is 
just  the  binomial  coefficient 


(i,-33)  I     ]   "  NFD!  (m-NFD)! 

\NFD/ 


Hence, 


(4.34)   IZ<u|vJu»>  =Y1     ■(  ](1-P)m-NFDPNFD  =  [(1-P)+P]m  =  1 

u  NFD=0  VNFD/ 

using  the  binomial  theorem.   Consequently,  V   is  stochastic, 
so  that 


(4.35)  V1  <=>A     (stochastic  case) 
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One  important  question  remains.   How  does  one*  move  a 

particle  according  to  the  conditional  probability  (V.  )   , ; 

1  uu" 

that  is,  given  a  particle  is  at  u',  how  does  one  sample  (V, )   , 
for  a  new  position  u?   The  following  routine  will  do  the  job: 

Independently,  and  with  probability  P,  flip  each 
of  the  m  spins  in  |u'^>  (=  \a  '  ,o'  .  .  .a'  >  ).   The 

resulting  ket  |u^>  is  drawn  from  (V,  )   ,. 
&      i  /  l'uu' 

The  proof  is  simple.   From  equation  (4.32), 


(4.36)  (Vl)uu,  =  pNFD(1_p)m-NFD   ^ 


but  the  right-hand  side  of  equation  (4.36)  is  just  the  probability 
of  moving  from  |u'^>   to  lu/1  by  the  above  routine. 

2 .  An  Improved  Walk 
The  walk  outlined  in  the  previous  section  proved  un- 
satisfactory.  Lattices  no  larger  than  4  x  4  x  «>  (three 
dimensions)  could  be  simulated  because  of  the  poor  stability 
and  accuracy  of  the  walk.   The  major  reasons  for  this  can  be 
seen.   The  fact  that  each  spin  in  |u'^>  was  flipped  independently 

means  that  most  of  its  correlation  with  both  other  spins  and 

B 
the  magnetic  field  was  supplied  by  the  fissioning  factor  —  , 

X 
after  the  move  was  completed.   Therefore,  if  a  particle  happened 

to  move  to  a  part  of  space  with  a  large  value  of  y     a. a..-,    and 

)>  a .  ,  it  produced  many  daughters.   Conversely,  if  it  moved  to 

an  unproductive  part  of  space,  it  produced  no  daughters.   From 
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the  point  of  view  of  importance  sampling,  this  is  just  what 
should  be  avoided. 

By  use  of  importance  sampling,  the  following  improved  walk 
was  devised.   Given  the  ket 

(4.37)  u'  =  |a^,a^,...a^> 

build  the  ket 


(H . 38)  u  =  I  a  a  . . .a 


> 


lw2"   m 


in  the  following  way: 

1.  Flip  the  first  spin  such  that  ^a-,^>  =  M,  where 
M  is  a  guess  of  the  system  magnetization. 

2.  Flip  the  next  spin,  to  the  value  o.t    conditional 

on  a!  and  the  result  °.    , ,  in  a  manner  to  be  explained. 
1  i-l'  ^ 

3.  Repeat  step  2  for  all  but  the  last  spin. 

4.  Flip  the  last  spin  to  the  value  aM  conditional  on 
aM,  and  the  results  aM_-,  and  a-,. 


It  can  be  seen  that  this  procedure  can  build  into  the  jump  not 
only  a  preference  for  the  spins  to  align  with  the  magnetic  field, 
but  also  a  preference  to  align  with  themselves.   We  now  investigate 
this  in  detail.   Given 

(4.  39)  <ulvluV>  =  <ulv2V3V1lu')> 

=    (V0)    (V_)     (V, )       , 
v    2'uv    3uv    1'uu'    , 
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being  careful  to  obey  the  commutation  relations  and  explicitly 
using  periodic  boundary  conditions,  define  the  following  "single" 
spin  operators 
(4.40)     <^u|V|u'>  = 

(Va  (V2V3)a  (Va  a»  '  ' '  (V2V35a  (Va  a'^Wa  a' 
d   al  d    3   °m  ram       J  °2   x  °2°2   J  °1    °1°1 


where 


(4.41)     (V1)aa,  =  (/a.|e-K+eKcx|a!>, 


1  i 


(4.42)  (v3)a.    =   exp[hai]    ,  ^^i-l^a      =   expEKc^a^]    , 


(4.43) 


a,   =   °M 


For  step  one  above,  multiply  and  divide  the  matrix  element  by 
P„   ',  defined  such  that 

<Vi 


(4.44) 


<°1> 


P    i    =  M  ; 
a1o1 


that  is, 


(1.15) 


■■•'¥A<Vo;<V» 


2  2 


1    P 


a-,a 


11 


alal 


For  step  2  do  the  following:   Define  the  (unnormalized)  probability 
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(4.46) 


Pr^la^V    E    'VVVa^Va^ 


C(a',a..)    is    its    normalization 


(4.47) 


cCoj,^)  =  YZ   ^3)^)   , 


Again  divide  and  multiply  the  matrix  element 


(4.48) 


C(o^a1)(V2V3)a2(V1)a2CT2]Pr{a2|a'a1} 


Pr{a^  \a'2,o-L} 


C(a^,a1) 


Vi 


SP 


alal 


On  simplifying  this  becomes: 


(4.49)  ^u|vlu/ 


[Pr^lojv"! 


C(a2a1) 


C(a2a1) 


^Vi 


alal 


alal 


where  the  two  terms  in  curly  brackets  become  the  flipping 

probabilities,  and  the  other  two  terms  become  the  developing 

fission  factor. 

Step  3  can  be  seen  to  be  a  generalization  of  step  2,  in 

which  the  flipping  probability  for  the  m'th  spin  becomes 

conditional  on  o1 ,o      , ,  and  a, . 

m'  m-1'       1 

By  use  of  this  revised  scheme,  the  efficiency  of  the 
computation  was  increased  by  a  factor  of  magnitude  10J  relative 
to  the  first  walk,  and  was  used  in  the  final  calculations. 
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D.  Estimating  the  Eigenvalue  and  Its  Logarithm 
The  random  walk  as  described  in  the  previous  sections  will 
provide  a  population  of  points  whose  distribution  is  T,  the 
"largest"  eigenvector  of  the  transfer  matrix.   What  is  desired 
is  the  largest  eigenvalue,  A,  which  is  the  partition  function 
per  layer,  and  the  free  energy  per  spin,  P(N*): 


k  T 
(4.50)         F(N*)  =  -  -jL-  In    A, 


N* 


where  N*  is  the  number  of  spins  per  layer. 

In  the  previous  sections,  it  had  been  assumed  that  X  was 
known.   This  not  being  the  case,  suppose  a  trial  eigenvalue,  A™, 
is  used.   In  Appendix  A  it  is  shown  that  the  mean  number  of 
daughters  per  point  is  unity  if  Arn  =  A.   Since  t —  is  a  constant 
multiplying  V|y^>  ,  this  statement  may  be  generalized  to 


/,,,--,  x      /  #  of  daughters  x>     ,  ,. 

(4.51)      <.  r2! y        =  A/Am  , 

\     points  /  T  ' 


or 


t  h    c^         x    ->  rx       r,  _  /  #  daughters 
(i<-52)         X    =    XTD  5      D  =  \ point  .' 


Therefore  a  Monte  Carlo  estimator  of  A  is: 


( ^.53)  AEST   AT  DEST 


where  D^^the  estimator  for.  the  mean  number  of  daughters  per 
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point,  can  be  scored  in  the  following  manner.   Let  P.  be  the 
number  of  points  in  the  i'th  iteration.   Assume  that  the  i'th 
to  the  (i+m)'th  iterations  are  to  be  scored;  then  an  estimator 
is  the  sample  mean 


(4.54) 


DEST 

l+m 

>—      PK 

K=i+1  A 

i+m-1 

>     PK' 

K'=i 

If  the  population  is  distributed  according  to  V,    i.e.,  fully 
converged,  this  estimator  is  unbiased.   However,  a  better 
estimator  (lower  variance)  would  be: 


(4.55)  AEST  =  XT  YZ   (Puu.)±/Q  , 


where  Q  is  the  total  number  of  points  that  moved  in  the  above 
block  of  iterations  (see  discussion  concerning  equation  (4.15)). 
This  was  the  estimator  used. 

Perform  the  above  procedure  for  several  different  blocks 
of  iterations,  obtaining  several  different  estimates,  and  form 
their  sample  mean,  IVU™.   The  result,  i.e.,  the  random  variable, 
MFC,„,  with  mean  D  and  unknown  variance,  will,  by  the  central 
limit  theorem,  be  assumed  to  have  a  normal  distribution. 

To  obtain  the  free  energy,  equation  (4.50),  an  estimate  of 
Jin  A.,  is  needed.   Now,  the  mean  of  the  random  variable,  Jin  MF_„, 
is  not  equal  to  in   D,  but  as  shown  in  Appendix  B,  the  effect  of 
the  bias  can  be  neglected. 
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Therefore,  the  following  estimator  for  the  free  energy 
was  used: 


(4.56) 


•FEST(N*> 
kBT 


w   [In    \T+ln   MEST] 
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V.   A  Technical  Description  of  the  Programs 

Because  the  programs  for  the  two-and  three-dimensional 
Ising  models  differ  only  in  the  meaning  of  a  layer  and  the 
geometry  of  the  interactions,  only  the  program  for  the  two- 
dimensional  model  will  be  described. 

Working  in  a  direct  product  (az)  representation  of  dimen- 
sionality 2  ,  basis  kets  (row  configurations)  were  represented 
by  computer  words,  each  of  the  first  m  bits  of  a  word  re- 
presenting the  value  (±1)  of  its  respective  spin.   For  example, 
for  a  f ive-spin-per-row  lattice,  a  possible  basis  ket  would  be 

(5-1)  |+4"H  +  >   , 

and  would  be  represented  by  the  computer  word 

(5-2)  KET  =  10100. . . 

Since  a  sixty-bit  computer  was  used,  any  value  of  m  less  than 
sixty  one  could  be  accommodated  without  further  bookkeeping. 
Routines  were  created  for  the  following  purposes: 

1.  To  create  a  random  ket  from  a  given  distribution. 

2.  To  flip,  with  known  probability,  the  j ' th  spin  in  a 
given  ket (0<j<m) . 

3.  To  find  the  number  of  up  spins  (a.=+l)  in  a  given  ket 

J 

4.  To  find  the  number  of  positive  bonds  (o-a.+-]=+l)  in  a 
given  ket. 
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5.  To  sense  the  direction  of  the  j'th  spin. 
To  create  an  initial  population  sampled  from  c|>  ,  the 
positions  (u)  of  approximately  two  thousand  points  were  drawn 
from  a  distribution  in  u  space,  in  which  each  spin  component 
had  the  same  independent  probability  of  being  up;  that  is 

(5.3)         (<|>0)u  =  p(a1)p(a2).  .  .p(am)  ,     a   =  ±1 


where 


(5.4)         p(o  =+1)  =  1  -  p(a  =-1)  ,    all  j 

J  J 


By  equation  (2.13),  P  should  be  chosen  close  to  unity  for  a 
region  of  temperature  and  field  in  which  the  magnetization  per 
spin  is  expected  to  be  close  to  one,  and  should  be  chosen  close 
to  one  half  in  a  region  in  which  the  magnetization  per  spin  is 
expected  to  be  close  to  zero. 

As  described  in  Section  IV.  A,  each  point  in  the  population 
was  then  independently  moved  from  its  position  |  u  '^>  to  its  new 
position  \u)   by  flipping,  with  known  probability,  each  spin  in 
the  ket  |u'^>  .   The  point  was  then  split  into  £  points  by 
sampling  a  K   from  g   ,(?).   This  was  done  in  the  following 
manner.   Suppose  that  when  the  point  was  moved  from  |u'^>  to  |u)>  , 
the  numerical  value  of  the  fissioning  part  of  the  matrix  was 


F   , .   Let  I  be  the  integral  part  and  F  be  the  fractional  part 

of  F   , :  i.e., 
uu"      ' 
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(5-5)  F   ,  =  I  +  F   ; 


with  probability  F   let  £  =  I  +  1; 
with  probability  1  -  F   let  £  =  I  . 
Then  ^  is  a  random  variable  with  mean 


(5-6)   C?>=  Fr(I+l)  +  (l-Fr)I  =  I  +  Pr  =  Fuu,  . 

Therefore  £  satisfies  equation  (4.15). 

Successive  Monte  Carlo  iterations  were  performed  until 
stability  was  reached.   During  this  convergence,  the  trial 
eigenvalue  was  adjusted  in  the  following  way: 


NA 
(5.7)        XT  -  —  AT  , 


where  N.  is  the  number  of  points  after  the  move  and  Ng  the 
number  before  the  move.   Convergence  to  V  was  assumed  completed 
after  it  was  noticed  that  Am  had  small  fluctuations  about  a 
relatively  stable  value.   Arrival  at  V      was  further  spot  checked  by 
performing  a  very  long  chain  of  iterations,  thereby  watching  the 
long  term  behavior  of  A  . 

After  convergence  was  completed,  Am  was  held  fixed  and  the 
iterations  were  continued.  Equation  (4.56)  was  used  as  a  Monte 
Carlo  estimator  of  the  free  energy. 

To  check  the  Monte  Carlo  results,  a  numerical   procedure 
was  developed  for  solving  the  square,  three-spi n-per-row  problem 
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(m=3)  •   The  resulting  eight-by-eight  transfer 
matrix  was  diagonalized  by  standard  numerical  methods.   Both 
above  and  below  the  critical  temperature  the  Monte  Carlo  result 
agreed  with  the  numerical  result  to  within  a  few  parts  in  ten 
thousand.   Also,  the  error  was  always  the  same  order  of  magnitude 
as  the  estimated  standard  deviation. 
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VI.   Results 

Using  the  previously  described  Monte  Carlo  procedure,  the 

zero-field  free  energy  was  found  for  simple  cubic  lattices  of 

12  v 
various  sizes.   Also,  the  exact  series  expansion  (to  v   )  of 

Wakefield  was  summed  at  the  same  temperature  to  give  a  good 

approximation  of  the  bulk  free  energy.   The  "temperature" 

chosen  was  Kp/l.83,  high  enough  that  the  series  in  v(tanh  K=.17) 

was  rapidly  converging.   Table  1  contains  the  numerical  results, 

wherin  the  interval  of  error  is  twice  the  estimated  standard 

deviation.   A  plot  of  free  energy  versus  lattice  size  is  shown 

in  Figure  1.   The  effect  of  surface  terms  (equations  2.48  and 

2.49)  can  be  seen  by  subtracting   them  from  the  Monte  Carlo 

results  (Fig. 2 ) . 

As  explained  in  Appendix  C,  the  magnetization  and  spin 
correlation  were  also  obtained.   Figure  3  shows  a  sample  magnetic 
isotherm  at  "temperature"  Kp/1.25,  slightly  above  the  critical 
point.   The  corresponding  numerical  data  are  contained  in  Table  3 

Defining  p(r)  as  the  spin  correlation  at  zero  field 


(6.1)  p(r)  =  <^i+r>  , 


in  which  both  spins  are  in  the  same  row,  both  the  nearest 
neighbor,  r  =  1,  and  the  next-nearest  neighbor,  r  =  2,  spin 
correlations  are  shown  in  Figure  4  for  various  temperatures  down 

to  the  critical  point.   Table  4  contains  the  corresponding 

2 
numerical  results.   A  plot  of  p  (1)  versus  p(2)  is  shown  in 

Figure  5. 
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TABLE  1 


Free  energies  of  the  simple  cubic  lattice 


N* 

-F(N*)/kDT 

D 

-Fr(N*)/kBT 

-VkBT 

9 

.717270  ±  1.9X10-5 

.715880 

.71582190^9 

16 

.715951  ±  l.lxlO-5 

.71580^ 

25 

.71586^  ±  4.7><10~6 

.715847 

49 

.715828  ±  5.^xl0~6 

.715828 

121 

.715822  ±  6.3X10-6 

.715822 

F(N*)  is  the  free  energy  per  spin  for  a  lattice  with  N* 

spins  per  layer. 
F  (N*)is  the  same  free  energy  with  all  surface  terms,  to 

order  m  +  4,  removed. 
F     is  the  bulk  free  energy. 


^7 


TABLE  2 


Magnetization  per  spin  for  cubic  lattice  with  121  spins  per  layer 


tanh(h) 

M  ( h ) 

.025 

.1955  ± 

.005 

.05 

.3295  ± 

.01 

.1 

.4898  ± 

.015 

.2 

.6456  ± 

.004 

.4 

.8156  ± 

.007 

K  =  Kc/1.25 

=  .17737 

TABLE  3 


Spin  correlation  for  a  cubic  lattice  with  121  spins  per  layer 


K 

P(D 

p(2) 

K'^1.83  = 

.12115 

.12979  ± 

.00098 

.01597  ± 

.0010 

K^l.50  = 

.14781 

.16284  ± 

.0012 

.02896  ± 

.0015 

Kc/1.25  = 

.17737 

.20676  ± 

.0031 

.05384  ± 

.0027 

Kc  = 

.22171 

.30557  ± 

.011 

.12830  ± 

.017 
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figure  2      The  effect  of  surface  terms 
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figure    3        Magnetization    Per   Spin 
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figure   4      Spin    Correlation 
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figure    5 
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VII.   Conclusions 

The  significance  of  the  numerical  results  will  now  be 
discussed.   The  exponential  form  of  the  free  energy  per  spin 
versus  lattice  size  (Figure  1)  shows  excellent  agreement  with 
the  result  obtained  from  high  temperature  graphical  analysis, 
equation  (2.42).   In  Figure  2  it  is  shown  that  the  low  order 
surface  terms  dominate  the  difference  in  free  energies  for  the 
lattices  with  m  equal  to  three  and  four.   However,  this  is  due 
only  to  the  particular  temperature  chosen;  the  surface  terms 
will  become  more  dominant,  though  smaller,  as  the  temperature  is 
raised.   The  sample  magnetization  isotherm  (Figure  3)  at  a 
temperature  above  but  near  the  critical  temperature  shows  a 
large  slope  (susceptibility)  near  h  =  0.  The  spin  correlations 
for  zero  field, 


p(r)  =  <a1a.+r>  , 


along  the  same  row  (or  column)  are  shown  in  Figure  4  for  r  equal 
to  one  and  two.   The  range  of  temperature  is  from  Kc/1.83  (high 
temperature),  down  to  the  critical  point,  Kc .   By  plotting  p  (1) 

versus  p(2),  Figure  5  suggests  exponential  behavior  which  is 

24 
expected  from  high  temperature  graphical  expansions. 

The  Monte  Carlo  method,  using  the  transfer  matrix  approach, 

has  been  shown  to  be  a  successful  way  of  solving  large  lattice 

problems.   At  present,  work  is  being  done  to  test  the  static 
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25 
scaling  hypothesis    for  finite  lattices,  and  to  include  in  the 

Hamiltonian  next  nearest  neighbor  interactions.   The  first  of 
these  tasks  will  prove  the  easiest,  graphing  various  magnetization 
isotherms  in  terms  of  reduced  variables.   The  second  task,  while 
straightforward,  requires  a  major  change  in  the  computer  program. 

This  thesis  will  conclude  with  a  brief  review.   Starting 
from  the  nearest  neighbor  Ising  Hamiltonian,  and  using  the 
canonical  ensemble,  the  partition  function  for  a  finite  lattice 
was  found  to  be  the  largest  eigenvalue  of  a  matrix  called  the 
transfer  matrix.   The  dimensionality  of  the  matrix  varied  ex- 
ponentially with  the  "size"  of  the  lattice.   Although  the  theory 
was  developed  in  terms  of  the  square  lattice  of  "size"  m  x  oo, 
the  actual  calculations  were  performed  suing  the  cubic  lattice 
of  "size"  m  x  m  x  <».   The  largest  eigenvalue  was  shown  to  be 
non-degenerate  and  the  corresponding  eigenvector  contained  only 
non-negative  elements.   The  high  temperature  power  series  ex- 
pansion of  the  partition  function  was  obtained  by  graphical 
techniques.   The  lowest  order  difference  between  the  free  energy 
of  the  bulk  lattice  (infinite  in  all  directions)  and  the  finite 
lattice  was  shown  to  consist  of  surface  terms  whose  corresponding 
graphs  circumscribed  the  lattice. 

'After  introducing  the  general  theory  of  Monte  Carlo,  a 
correspondence  was  established  between  the  eigenfunction  belong- 
ing to  the  largest  eigenvalue  of  the  transfer  matrix  and  the 
steady  state  collision  density  of  a  class  of  random  walks.   After 
choosing  a  particular  random  walk,  a  method  of  obtaining  the 
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largest  eigenvalue  was  found,  which  consisted  of  statistically 
sampling  the  walk.   Knowledge  of  the  eigenvalue  produced  the 
partition  function  directly,  and  after  an  analysis  of  the  bias, 
produced  the  free  energy.   By  removing  the  surface  terms  from 
the  free  energy,  and  comparing  this  reduced  free  energy  to  the 
bulk  free  energy,  the  effect  of  the  finite  size  of  the  lattice 
was  shown.   A  method  of  obtaining  numerical  derivatives  of  the 
partition  function  was  found,  which  consisted  of  associating  a 
statistical  weight  to  each  point  in  the  walk.   Using  this  method 
both  magnetization  isotherms  and  spin  correlations  were  found. 
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Appendix  A.   The  Branching  Process 

Proof  that  the  splitting  technique,  as  defined  below, 
produces  a  random  walk  whose  density  is  <j> ,  where 


(Al)  4>  =  ?£  <j> 


X    is  the  eigenvalue  of  BA , 

A  is  the'  stochastic, 

B  is  positive, 
and  that  the  mean  number  of  daughters  per  point  is  unity.   Rules 
of  the  walk: 

1.  Given  a  point  at  u'  drawn  from  p.d.f.  <f>  ,,  jump  it 

to  u  by  means  of  the  conditional  probability  A   , . 
j  i  a      uur 

2.  From  a  density  g   ,  (£)j  where 


<^C>     =    B   ,/X 
X  /     uu' 


sample  a  £  . 
3.  Split  the  one  parent  point  into  E,     points,  each  of 
which  remains  at  u.   For  example,  if  5  =  3>  the 
single  point  at  u  becomes  three  points  at  u. 
In  any  random  process,  the  mean  number  of  outcomes  of  an 
event  £  equals  the  number  of  trials  times  the  probability  of 
the  outcome  &.   For  example,  the  mean  number  of  heads  in  N  flips 
of  a  coin  is  N/2.   Using  this  principle,  consider  N  points  in 
a  population,  the  position  u  of  each  drawn  from  density  cj> 


u 
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Move  each  point  independently  according  to  the  above  rules.   Let 

C(u')  be  the  mean  number  at  u'  before  the  move. 

J(u,u')  be  the  mean  fraction  of  these  moving  to  u. 

D(UjU?)  be  the  mean  number  of  daughters  created  per 

arrival  at  u  from  u'. 

M(u,u')  be  the  mean  number  arriving  at  u  from  u'. 
Then, 

(A2)  M(u,u*)  -  D(u,u' )J(u,u' )C(u'  )  . 

Let  M(u)  be  the  mean  number  arriving  at  u  from  any  u'. 

(A3)  M(u)  =  YL"   M(usu')  =  JZ  '   D(u)J(u,u')C(u')  . 

u1  uf 

But,  by  definition, 


l^)  D(u,u»)  =   Buu,/X  ; 


also 

(A5)  C(u»)  =  N(f>u, 

N  being  the  population  size  before  the  jump,  and 


J(u,u' )  =  A   , 
'       uu' 
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Therefore,    comparing   equation    (A3)    with   equation    (Al), 


M(u)    =   N<[>u      , 


proving  that  the  population  stochastically  duplicates  itself 
after  the  moves.   Since  by  definition  §   is  normalized, 


(A8)  T~  d>u  =  i  , 


the  mean  number  of  points  (N)  after  the  jump  equals  the  number 
of  points  before  the  jump  (N),  since 


(A9)  N  =  YH   M(u)  =  EZ  N({>   =  N 

u  u 


Therefore,  the  use  of  the  correct  eigenvalue  \    stabilizes  the 
population.   This  fact  is  used  to  provide  a  Monte  Carlo  estimator 
of  X  . 
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Appendix  B.   The  Mean  and  Variance  of  a  Lognormal  Estimator. 
Let  X„  be  a  random  variable  with  normal  density: 

-(X  -X)2 
(Bl)  P(XF)  =  exp[  ~—  ]  , 

^4  2QxE 

E 
2 
where  Q  ,  is  the  variance. 

A 

E 
Expand  Jin  Xp  in  a  Taylor  series  about  X: 

(B'2)  In   XE  =  Jin  X  +  yCXg-X) ^(X£-X)2+  -^(Xg-X)3-  ORXg-X)1*] 

2X  3X 


Take  the  mean  of  Jin  X_  over  the  normal  distribution 

E 


(B3)  M(£n    XE)    =  Tin    Xg    P(XE)dXE 


Substituting  the   expansion   for   Jin   XE, 


(B4)         M(Jln   XE)    =    Jin   X ^"     /  (Xg-X)  2P(X£)dXE 


The  second  and  fourth  terms  in  the  expansion  give  zero  contribution 
Therefore , 

(B5)  H(nn  Xc)  =tnl- 


E'    "•      2X2 


proving  that  &n  X„  is  a  biased  estimator  of  £,n  X,  with  bias  B. 
(B6)  B  = 


2 

2X^ 
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to  fourth  order. 


The  variance  of  In   X^  may  be  similarly  obtained.   Expand 


£n2XT 


,(B7)   £n2XT7  =  £n2X  +  2£^  X(X„-X)  + 


(  Y     -Y  V 

2_    2&n  X  |   v  E  ' 

2       2 

X^      X^ 


Again,  the  first  and  third  order  terms  give  zero  contribution: 


(B8) 


QX     Qx  ln 

M(£n2XTT)  =  £n2X  +  — § ^-5- 

**  X^        X^ 


Now, 


(B9)  V(£n  XE)  =  M(£n2XE)  -  [M(£n  Xg)]2 


where  V  is  the  variance. 


QX 
(BIO)   VUn  XE)  =  £n2X 1 


2     rt2 


Qv  £n  X 
XE 


A 


2 


Qx  in   X 
£n2X ^~ 


/^ 


'Q 


X 


2 


XE 

7. 


so  that 


(BID 


Q 


XT 


V(£n  XE) 


X" 


or 


(B12) 


Q 


in   XT 


X 
"X 


E 
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For  any  estimator,  X  ,  worthy  of  the  name,  and  in  particular 
for  the  results  in  this  thesis, 

QX 
(B13)  -jp  <  1  ; 

but,  from  equation  (B6), 


(B14)  B  "  laL  XF  <K    a£n  Ap 

Ei  Ei 


for  any  good  estimator.   Therefore,  the  bias  of  the  estimator, 

£n  X_,  is  so  much  smaller  than  the  measure  of  statistical  error, 
E'  ' 

a„   „.  that  it  can  be  neglected. 

Jen  E5  & 
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Appendix  C.   Magnetization  and  Spin  Correlation. 

To  obtain  both  the  magnetization  and'  the  spin  correlation, 
logarithmic  derivatives  of  the  free  energy  are  needed.   From 
equation  (1.15),  the  magnetization  per  spin  is 


/m\  m/i~\         1  d    In    Z  1  3  In    X 

<C1)         M(h)  =  N—Ih-  =E-^h—  ' 


where,  as  previously,  the  square  lattice  is  being  considered 
for  convenience  of  notation. 

The  spin  correlation  may  be  defined  as 


(C2)  P(r)  =  (^^i+r/ 


for  two  spins  in  the  same  layer,  separated  by  a  "distance"  r. 
This  function  can  be  put  in  the  form  of  a  logarithmic  derivative 
of  the  free  energy  by  means  of  the  following  technique.   In  each 
layer  of  the  lattice,  create  a  coupling  between  spin  i  and  spin  i+r. 
That  is,  add  to  the  Hamiltonian  the  term 


N  _n 
(C3) 


"kBTL  C  °i,J*i+r,J 


where  i  is  the  intralayer  index  and  j  is  the  interlayer  index, 
The  new  Hamiltonian  becomes 


n 


(do  4]'  ="H-I  -  ^Cvi+w 


6h 


then. 


(C5) 


and 


(C6) 


exp 


U-,  }Up  . 


n 


i4 


k  T 
KB 


exp 


j=l 


La .  . a .  . 


3  g,n  2' 
9L 


n 


L=0 


a .  .a. , 


where  the  average  is  for  the  standard  Ising  lattice 
the  translational  symmetry  of  the  lattice, 


Using 


n 


(C7) 


and 


(C8) 


"5   o.  .a .  . 
fry   i,J  i+r,J 


n /  a. ,a 


L+r  , 


P(r)  = 


1  9  in   Z' 

n    3L 


L=0 


producing  the  result 


(C9) 


/  >,  _  9  £n  A 
P(r) 3L— 


L=0 


Since  both  M(h)  and  p(r)  are  logarithmic  derivatives  of  the  free 
energy,  the  Monte  Carlo  methods  of  solution  will  be  developed  in 
terms  of  the  magnetization. 

Using  the  finite  difference  approximation 


(CIO) 


M(h)  : 


1  r£n  x(h+Ah)-£n  \(h)   -i 
Ah  J 


rn 
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and  assuming  'cnat  Dutn  in   Ath+A)  and  la  A(h)  are  exactly  known, 

the  error  involved  is  negligible,  because  Jin  X   is  a  smooth 

v 
function  of  h  away  from  the  critical  point,  and  Ah  will  be  of 

-7 
order  10 

In  standard  numerical  procedures,  use  of  the  finite 
difference  approximation  should  be  avoided  if  possible  because 
any  small  error  in  either  in   A(h+A)  or  X(h)  produces  a  large 
errbr  in  the  finite  difference  derivative.   In  fact,  this  error 
gets  larger  as  Ah  gets  smaller.   Fortunately,  Monte  Carlo  avoids 
this  problem.   If  the  Monte  Carlo  estimate  of  in   X(h+Ah)  is 
highly  correlated  with  the  estimate  of  in   A(h),  any  error  in 
one  will,  with  high  probability,  occur  in  the  other  and  will 
cancel  in  the  difference,  in   A(h+A)  -  Jin  X(h). 

This  correlation  was  produced  by  the  use  of  statistical 
weights,  which  will  now  be  described.   After  stability  has 
been  reached,  associate  with  each  point  (i)  in  the  population 
a  number  w . .   Set  each  w.  equal  to  1.   Iterate  the  population 


in  the  standard  fashion, 

3       V  (h)V  V   -, 
(Cll)         Y(h)  =  -* £-i-  Y(h) 


where  the  dependence  on  h  is  shown  explicitly.  After  each  point 
moves  from  its  old  position,   u1   ,  to  its  new  position,   u 
multiply  the  point's  weight  by 


V  (h+Ah)   . 
(C12)         <u|   3        |u^> 
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i.e., 


(C13) 


w.  =  w. 

1 


V  (h+Ah) 


i  C1  ivTTh) 


■y 


After  fissioning  the  point,  associate  this  new  weight  with 
each  of  its  daughters. 

After  the  iteration,  the  particle  density  is,  as  usual, 
V(h).      However,  the  weight  density  is  Y(h+Ah),  as  will  now  be 
shown.   Let  the  mean  arrival  of  weight  at  u  be  %.      Then, 


(C14) 


^u  =  I 


/V3(h+Ah) 


u' 


V 


V-(h) 
3    /uu 


(V  (h))   (v2)   (V  ) 

'   uu    uu uu  '  .j 


u 


In  words : 


The  mean  arrival  of  weight  at  u,  %   ,  equals  the 
mean  weight  at  u',  %. ■  ,,  times  the  mean  fraction 
moving  from  u'  to  u,  times  the  mean  number  of 
daughters  per  point,  times  the  number 


(C15) 


V3(h+Ah) 
V3(h) 


uu 


summed  over  all  u1 
Simplifying, 


(VJ   (Vn) 


(C16) 


Xu  =  ZZ  (V.(h+Ah)) 
u  t    j 


uu    uu 


uu 


X 


But  this  is  the  equation  satisfied  by   (h+Ah).   Therefore, 


(C17) 


3       => 

*%  =  Y(h+Ah) 
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The  result  is  that  the  bare  points  are  distributed  according  to 
tCh),  and  the  same  points  with  associated  weights  are  distributed 
according  to  YCh+Ah). 


The  difference  in  equation  (CIO)  which  must  be  estimated  is 


(C18) 


R  =  In    X(h+Ah)  -  In    X(h) 


To  estimate  in   A(h),  equation  (4.56)  was  used.   The  estimator  of 
In    A(h+Ah)  is  just  an  extension  of  equation  (4.5*0  to  the  case 
of  weighted  points: 


(C19) 


in   A(h+Ah)EST  =  £n  A?  +  in 


(WA 


W 


B 


EST 


where 

W„  =  sum  of  the  weights  after  the  iteration, 
W„  =  sum  of  the  weights  before  the  iteration. 

Sample  magnetization  and  correlation  curves  are  contained  in 

Section  VI. 
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